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ABSTRACT 



Eclipsing binaries are extremely attractive objects be- 
cause absolute physical parameters (masses, luminosi- 
ties, radii) of both components may be determined from 
observations. Since most efforts to extract these pa- 
rameters were based on dedicated observing programs, 
existing modeling code is based on interactivity. Gaia 
will make a revolutionary advance in shear number of 
observed eclipsing binaries and new methods for auto- 
matic handling must be introduced and thoroughly tested. 
This paper focuses on Nelder & Mead's downhill simplex 
method applied to a synthetically created test binary as it 
will be observed by Gaia. 

Key words: binaries: eclipsing; methods: numerical, sta- 
tistical; missions: Gaia. 



1. INTRODUCTION 



Many assessments have already been done for Gaia ex- 
pected harvest of eclipsing binaries (EB) to V<15, where 
both photometric and RV observations will have been 
available (see Munari et al. 2001, Zwitter et al. 2003, 
Marrese et al. 2004 and others). Out of 50 million ob- 
served stars, roughly 100000 will be double-lined eclips- 
ing binaries. However, based on experience from Hip- 
parcos, out of 1 billion stars observed to V<20, there 
will be ^2 million (Eyer, these proceedings) eclipsing bi- 
naries without spectroscopic observations, but with quite 
decent photometric accuracy. This study presents a new 
approach developed for automatic reduction of observed 
data along with an estimate of how much we may expect 
to obtain from them. 



2. THE METHOD 



Obtaining physical parameters from observations is an 
inverse problem solved numerically by a modeling pro- 
gram. Affirmed and most widely used is the WD code 
(Wilson & Devi nnevll 19711) . which features Differential 
Corrections method (DC) powered by the Method of 



Multiple Subsets (MMS) ( Wi lsonll993l) . DC method has 
already been applied successfully to automatic parameter 
extraction (e.g. Wyithe & Wilson 2001, Wyithe & Wilson 
2002, Prsa 2003 and others), but its original philosophy 
is based on interactive monitoring of each convergence 
step. The algorithm is very fast and works well if the 
discrepancy between the observed and computed curves 
is relatively small, but it tends to diverge or give physi- 
cally unmeaningful results if the discrepancy is large. As 
part of an effort to create a reliable and powerful package 
for EB analysis, a complementing minimization scheme 
is proposed. 



2.1. Nelder & Mead's Downhill Simplex 

There are two main deficiences of the DC method that are 
especially striking. 1) Once a DC method converges to a 
minimum, there is no way of telling whether that mini- 
mum is local or global; even if it is local, the method is 
stuck and cannot escape. 2) The main source of diver- 
gence and the loss of accuracy in the DC algorithm is the 
computation of numerical derivatives. 

To circumvent these two problems, Nelder & M ead's 
downhill Simplex 1 method (Nelder & Mead 1965), here- 
after NMS, is implemented. Since NMS doesn't compute 
derivatives but relies only on function evaluations, it is a 
promising candidate for our purpose. Basic form of the 
NMS metho d along with WD implem entation was first 
proposed by Kallrath & Linnell ( 1987). We take a step 
further and adapt the method specifically to EBs. 

NMS method acts in n-dimensional parameter hyper- 
space. It constructs n vectors p^ from the vector of initial 
parameter values x and the vector of step-sizes s as fol- 
lows: 

Pi = (x ,X 1 ,...,X i + 8i,...,X n ) (1) 

These vectors form (n + 1) vertices of an n-dimensional 
simplex. During each iteration step, the algorithm tries 
to improve parameter vectors pi by modifying the ver- 
tex with the highest function value by simple geometrical 

1 Nelder & Mead's downhill simplex is not to be confused with linear 
or non-linear programming algorithms, which are also referred to as 
Simplex methods. 
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transformations: reflection, reflection follo wed by expan- 
sion, contraction and multiple contraction (iGalassi et al.l 
2003). Using these transformations, the simplex moves 
through parameter space towards the deepest minimum, 
where it contracts itself. 



PHOEBE (http://phoebe.fiz.uni-lj.si) is a software pack- 
age built on top of the WD code that extends its basic 
functionali ty to encompass, am ong other extensions sum- 
marized in lPrsa & ZwitteJ ll2005l) . the NMS method. It is 
especially suited for EBs: powered by heuristic scans, pa- 
rameter kicking and conditional constraining, the method 
is able to efficiently escape from local minima. 



2.2. Heuristic Scan 

NMS is a robust method that always converges, but it 
can converge to a local minimum, particularly since pa- 
rameter hyperspace in vicinity of the global minimum is 
typically very flat, with lots of local minima. In adition, 
global minimum may be shadowed by data noise and de- 
generacy. 

Heuristic scan is a method by which a minimization algo- 
rithm selects a set of starting points in parameter hyper- 
space and starts the minimization from each such point. 
It then sorts all solutions by the cost function (the x 2 , 
for example) and calculates parameter histograms and 
convergence tracers for given hyperspace cross-sections 
(specific examples are given in Section 0. The way of 
how the algorithm selects starting points is determined 
by the user: the points may be gridded, stochastically 
dispersed, distributed according to some probability dis- 
tribution function (PDF) etc. The basic idea of heuristic 
scan is to obtain decent statistics of adjusted parameter 
values from which a fair and realistic error estimate may 
be given. 



where index i runs over M measurements within a sin- 
gle data-set and index k runs over N data-sets (different 
photometric curves). Since the variance is given by: 

4 = 1 Xfei -X) 2 : ( 3 ) 

i 

we may readily express xt as: 



x 2 = (N k -l)^ 



(4) 



and the overall % 2 value as 

X 



If <7fc are fair and all data-sets contain approximately the 
same number of observations, then the ratio Sk/uk is of 
the order unity and \ 2 of the order NM. This we use for 
parametrizing x 2 values: 

NM, A := (x 2 /xl) ■ quantization. (6) 



Parameter kicking is a way of knocking the obtained 
parameter-set out of the minimum: using the Gaussian 
PDF, the method randomly picks an offset for each pa- 
rameter. The strength of the kick is determined by the 
Gaussian dispersion <7kick, which depends on A: if the 
value is high, then the kick should be strong, but if it is 
low, i.e. around A ~ 1, then only subtle perturbations 
should be allowed. Experience shows that a simple ex- 
pression such as: 



Ckick 



0.5A 

Too 



(7) 



works very reliably. This causes Ckick to assume a value 
of 0.5 for lOcr offsets and 0.005 for ler offsets, being lin- 
ear in between. Note that this Ckick is relative, i.e. given 
by: 

(8) 



abs 
"kick 



T rr rel 
"kick! 



2.3. Parameter Kicking 



where x is the value of the given parameter. 



Another possible approach to detect and escape from lo- 
cal minima is to use some stochastic method like Simu- 
lated Annealing (SA). However, such methods are notori- 
ously slow. Since the EB hyperspace is typically very flat, 
stochastic methods would be practical only in the vicin- 
ity of the global minimum. Thus instead of full-featured 
S A scan, a simple new procedure has been developed that 
achieves the same effect as stochastic methods, but in sig- 
nificantly shorter time. 

The idea is as follows: whenever NMS reaches a min- 
imum within a given accuracy, the algorithm runs a 
globality assessment on that minimum. If we presume 
that standard deviations a k of observations are estimated 
properly and that they apply to all data points regardless 
of phase or nightly variations, we may use them for x 2 
weighting: 

M M 

*l =y2w k (xi -x) 2 = —^2(xi ~x) 2 , (2) 

l — l K 4 — 1 



2.4. Conditional Constraining 

Having purely photometric (LC) observations, it is im- 
possible to determine absolute physical parameters of the 
observed binary. However, if the distance to EB is mea- 
sured independently or if additional assumptions about 
the EB are set, even purely photometric observations can 
yield absolute values of physical parameters. If addi- 
tional constraints are imposed on the model from the out- 
side, the model is referred to as conditionally constrained 
(CC'd). 

Two CCs are immediately evident: 1) astrometric mea- 
surements on-board Gaia may be used to measure the dis- 
tance with t he accura c y of ~ 1 1 /Ltas at V=15 to -~ 165 /ias 
at V=20 fES A-SCfl <2000l) 41: 2) since a substantial 
number of stars are main sequence stars, M-L, L-T 
and T-R relations m ay be adopted as constraints. See 
iPrsa & Zwitterl d20 05) for details on CC implementation 
in PHOEBE. 
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3. SIMULATION 



To test the suitability of NMS for EBs, we built a 
partially-eclipsing synthetic main-sequence F8 V-Gl V 
binary using PHOEBE. Table ^ lists a U of its princi- 
pal parameters. The simulation presented in this paper 
is based exclusively on photometric data: two sets of 
observations corresponding to Johnson B and V filters 
are created, each with 82 points with Poissonian scatter 
Cobs = 0. 02mag at V = 20, values typical to expect 
from Gaia (ESA-SCI f200jj)4). 

Simulation flow is the following: all physical parameters 
set for adjustment (a, i, T\, T%, VL\, fi2> L\ and L2) were 
displaced by ^50%. The obtained set was used as ini- 
tial guess for the NMS. In the first part of the simulation 
the method converges to a solution using only heuristic 
scan and parameter kicking, which yields relative values 
of parameters. In the second part the simulation was addi- 
tionally CC'd by the main-sequence constraint to obtain 
absolute values of parameters. 
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Figure 1. \-histograms for 3 consecutive parameter 
kicks. Top-left figure demonstrates how numerous are lo- 
cal minima and how difficult it is for NMS to circumvent 
them. Other three figures show significant improvement 
by using parameter kicking. Histograms consist of 20 
bins (single bin width is 0.01). The last bin encompasses 
all higher values of A. Labels (a) through (d) are used 
consistently throughout the paper. 



Table 1. Principal parameters of the simulated binary. 



Parameter [units] 


F8 V 


Binary 


Gl V 


6 


P [days] 




1.000 




5 


a [R Q ] 




5.524 




4 


q = m 2 /mi 




0.831 




cf 


i [°] 




85.000 




7 


Teff [K] 


6200 




5860 




L [L ] 


2.100 




1.100 


6 


M [M ] 


1.236 




1.028 




R [Rq] 


1.259 




1.020 


5 


q [-] (a) 


5.244 




5.599 



O) Unitless e ffective potentials defined by 
IWilsonlff?^ . 




4. RESULTS 

Heuristic scan over all adjustable parameters has been 
performed to obtain accurate convergence statistics: over 
2 000 uniformly distributed starting points in parameter 
hyperspace were used during simulation. We present the 
results of overall minimization step-by-step. 

a) Convergence assessment. Depending on the bumpi- 
ness of the hyperspace, heuristic scan will generally yield 
different solutions from different starting points; it is our 
hope that only few of these solutions will account for 
most scans. To evaluate their quality, the globa lity as- 
sessment mechanism introduced in section 12.31 is used 
to sort solutions by the depth of the reached minimum. 
Tests show that the NMS method itself is all-too-often 
stuck in local minima and only ^15% of all runs end 
up within one percent of ideal A (A=l in case of Fig.Q. 
However, tests also show that parameter kicking signif- 
icantly improves this percentage (^50% after the first 



Figure 2. Convergence tracers for 3 consecutive param- 
eter kicks. These plots trace convergence steps within 2D 
cross-sections of the hyperspace, revealing areas of min- 
ima and degeneracy. Cross-hairs mark the location of the 
global minimum. 



kick, ^63% after the second and ^75% after the third 
kick); even more, parameter kicking also enhances con- 
vergence speed after each kick. 

Although the values of A may seem promising, they don't 
necessarily guarantee that the corresponding solution is 
optimal. Rather, additional assessments should be done. 
Fig. |3 shows convergence tracers: 2D cross-sections 
of parameter hyperspace tracing convergence from each 
starting point of heuristic scan to the corresponding so- 
lution. Such tracers clearly show areas of minima and 
degeneracy. Since the location of the global minimum is 
not known in real life, extra care should be taken to never 
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Figure 3. Histogram of the inclination for 3 consecutive 
parameter kicks. Histogram consists of 20 bins with 0.5° 
each. The solution is symmetric to i = 90°, but we adopt 
i < 90° by convention. 



blindly trust the statistics of such a degenerate problem. 

b) Statistics of obtained parameters. 

The usual practice in literature, when listing obtained pa- 
rameters from the model, is to give their values with for- 
mal errors, i.e. standard deviations reported by the used 
numerical method. These errors are often too optimistic, 
since degeneracy and noise noticeably contribute to the 
overall error. NMS powered by heuristic scan and pa- 
rameter kicking has the advantage of obtaining parame- 
ter errors statistically, independent of the method itself. 
Fig. [3] shows an example of obtained histogram for in- 
clination i. It is evident that for the NMS without pa- 
rameter kicking (and similarly for any other numerical 
method that cannot escape from local minima) any eclips- 
ing system is a tie; for a fully minimized solution (bottom 
right plot on Hig. 0» the error is simply standard deviation 
of the Gaussian being fitted over the histogram, yielding 
roughly 0.5°. 

c) Conditional constraining. 

Inclination is the only intrinsic parameter that may be ob- 
tained in absolute sense from photometric observations. 
Using CC this deficiency is removed: any particular CC 
adds one or more implicit parameter ties into the sys- 
tem. It basically introduces an intersection plane with the 
otherwise degenerate part of the hyperspace, thus elim- 
inating degeneracy. Fig. |4] demonstrates how a main- 
sequence constraint breaks the degeneracy for gravity po- 
tentials Oi 2- 



5. DISCUSSION 



The idea behind the NMS implementation is not to re- 
place the DC method but to complement it. DC is created 
for interactive usage and converges in discrete steps that 
need monitoring. NMS on the other hand aims to auto- 
mate this process so that intermediate monitoring is no 
longer necessary, which is a key goal for Gaia. DC is one 
of the fastest methods (WD's DC in particular, since it is 
optimized for EBs), but may easily diverge. On expense 
of speed, NMS is one of the most robust algorithms for 




Figure 4. Main-sequence constrained convergence trac- 
ers for 3 consecutive parameter kicks. Comparing this 
result to Fig.^\clearly shows that the intersection of both 
areas indeed gives the right solution. This is expected, 
since our synthetic binary is in fact a main-sequence bi- 
nary. 

solving non-linear problems and never diverges. Finally, 
both DC and NMS methods suffer from degeneracy and 
may be stuck in local minima. To overcome this, DC is 
complemented by the MMS and NMS is complemented 
by heuristic scan and parameter kicking. These differ- 
ences in intent make a combination of both methods a 
powerful engine for solving the inverse problem. 
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